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SPECIFICATION 



CROSS-REFERENCE TO RELATED APPLICATION 

This application claims the benefit of Provisional 
Application No. 60/116,362 filed January 19, 1999, 
incorporated herein by reference, 

BACKGROUND OF THE INVENTION 

This invention is concerned with signal processing of 
transients, and more particularly, with multiexponential 
signal processing. 

In the recent REVIEW article entitled "Exponential 
analysis in physical phenomena". Review of Scientific 
Instruments, Vol. 70, No. 2, Feb, 1999, pages 1233-1257 
(incorporated herein by reference) , authors Andrei A. 
Istratov and Oleg F. Vyvenko consider various aspects of 
multiexponential analysis. It is apparent from the text of 
the REVIEW article, and from the over 300 literature 
citations, that exponential analysis in physical phenomena 
is a subject of considerable interest in many scientific and 
technological disciplines, including, inter alia ^ solid 



state physics, medicine, biology and biophysics, geophysics, 
optics, engineering, chemistry and electro-chemistry. 

As pointed out in the REVIEW article, many physical 
phenomena are described by first-order differential 
equations whose solution is an exponential decay. The 
amplitude and time constant of the exponential decay carry 
information about the nature of the phenomenon being 
studied. As commonly happens, a number of exponential 
processes take place simultaneously, and equipment employed 
in analyzing such exponential processes (multiexponential 
decays) yields a signal which is a sxm of a plurality of 
exponential components . 

Obtaining useful information from multiexponential 
decays is not a simple task. Many methods of 
multiexponential analysis and many algorithms for this 
purpose have been proposed, but all of them have 
limitations . 

U.S. Patent No. 5,517,115 issued May 14, 1996 to 
Manfred G. Prammer (incorporated herein by reference) 
proposes a method and apparatus for efficient processing of 
nuclear magnetic resonance (NMR) echo trains in well 
logging. A priori information about the nature of the 
expected signals is used in an attempt to obtain an 
approximation of a model using a set of pre- selected basis 
functions. A singular value decomposition (SVD) is applied 
to a matrix incorporating information about the basis 
functions, and is stored off-line in a memory. During 



actual measurement, the apparatus estimates a parameter 
related to the signal-to-noise ratio (SNR) of the received 
NMR echo trains and uses it to determine a signal 
approximation model in conjunction with the SVD of the basis 
> function matrix. This approximation is used to determine, 
in real time, attributes of the earth formation being 
investigated. 

Techniques such as those described in the Prammer 
patent rely heavily upon a priori information about the 
0 nature of the expected signals. In attempting to obtain a 

reliable estimate of an unknovm model, i.e., parameters of a 
subject being investigated, such techniques rely on one or 
more numerical algorithms for multiexponential analysis 
which use so-called "fitting routines" in an attempt to fit 
.5 data to the model. 

Underlying the present invention is the discovery that 
substantially improved results can be attained without using 
fitting routines, and, instead, by emphasizing better linear 
resolution. As pointed in the REVIEW article, fitting 
2 0 routines work well only if the hypothesis of the n\jmber of 
multiexponential components is correct and an initial 
approximization is close to the true solution. The present 
invention, which does not use fitting routines, avoids the 
problems that are inherent in the use of such routines . 
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BRIEF DESCRIPTION OF THE INVENTION 

An object of the present invention is to provide 
improved methods and apparatus for the analysis of 
transients and for obtaining useful information therefrom, 
5 More particularly, an object of the invention is to 

provide improved multiexponential signal processing. 

One aspect of the invention involves a computer- 
readable medixna containing a set of coefficients that define 
a transform operator such as a matrix, 
0 Another aspect of the invention involves a method of 

calculating a transform operator utilizing a plurality of 
resolution functions , 

Another aspect of the invention involves a method of 
multiexponential signal processing in which multiexponential 
15 signals are sampled, and the above-mentioned transform 
operator is applied to the sampled signals. 

Another aspect of the invention involves an apparatus 
for multiexponential signal processing that comprises a 
signal processor including the above-mentioned transform 
2 0 operator* 

The present invention starts with the construction of 
appropriate transform operators. Once appropriate transform 
operators have been constructed, they are incorporated in 
signal processors of analytical instrumentation for 
25 processing data. Generally, such instrumentation includes a 
computer, as is well known in multiexponential analysis, 
Multiexponential data signals are sensed or detected by 

4 



conventional equipment and are input to the transform 
operator of the computer for signal processing • The signals 
may be applied in real time or they may be read out from a 
suitable storage medixim. 

5 Typically, the signals are applied to the transform 

operator in digital form after conventional sampling and 
analog- to-digital conversion. For example, digital samples 
of multiexponential decays may be obtained at equally- spaced 
instants in time, beginning at or just afte.r the start of a 

lO multiexponential signal* 

The invention makes use of linear resolution to obtain 
a better estimate of an unknown model. The present invention 
provides the very desirable property of optimal linear 
resolution of the unknovm model when plotted against the log 

15 of the time constant of the decay curves. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

The invention will be further described in conjunction 
with the accompanying drawings, which illustrate preferred 
(best mode) embodiments, and wherein: 
5 Fig. 1 is a diagram showing a matrix with coefficients 

• . - for transforming data d^ • . . dj, to produce 

parameters m^ . . . m^ of an estimate of an unknown model; 

Fig. 2 is a flow chart showing a method for calculating 
the coefficients of a row of a transform matrix in 
0 accordance with the invention; 

Fig. 3 is a table showing, inter alia , coefficients for 
three rows of a matrix, for r's of interest, namely 1, 10, 
100; 

Fig. 4 is a diagram showing data functions, resolution 
L5 functions, noise response, and point spread functions for a 
matrix constructed in accordance with the invention where 
the noise gain (NG) is 1.0000; 

Fig. 5 is a similar diagram for another matrix 
constructed in accordance with the invention for a noise 
20 gain of 3.1623; 

Fig. 6 is a diagram showing resolution func-tions for 
four matrices constructed in accordance with the invention 
with different noise gains; 

Fig. 7 is a diagram showing point spread functions 
25 associated with four matrices constructed in accordance with 
the invention for different noise gains; 
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Fig. 8 is a block diagram showing apparatus for 
processing data in accordance with the invention; 

Fig. 9 is a diagram showing decay curves of an MRI; and 
Fig. 10 is a diagram showing MRI relaxation 
distributions for four matrices constructed in accordance 
with the invention with different noise gains. 

DETAILED DESCRIPTION OF THE INVENTION 

One of the principal objectives of the present 
invention is to provide signal processing of transients, 
such as multiexponential decays, producing outputs that are 
better estimates of an unknown model to be investigated. 
Such estimates permit better interpretation of data, so that 
a user (researcher, physician, scientist, or engineer, for 
example) can obtain more accurate information as to the 
nature of the unknown model. Considered from one point of 
view, the invention may be looked upon as a better digital 
lens that provides improved resolution, just as a better 
optical lens provides improved resolution. 

A first step in achieving the objectives of the 
invention is to construct an improved transform operator, 
conveniently in the form of a matrix. The manner in which a 
transform matrix may be constructed, pursuant to the 
invention, will be described in detail later. The actual 
transform operator will depend upon its intended 
application, for example, medical imaging or well logging. 
For any application, several different transform operators 



may be constructed, to provide a user with greater 
flexibility. 

The present invention requires an understanding of what 
is referred to in the art as estimating a solution of a 
5 linear inverse problem. The linear inverse problem is one 

of communicating to an interpreter what is known and what is 
not known about an unknown model- The almost universal 
practice in the prior art for estimating the solution of a 
linear inverse problem is to calculate one or more of an 
.0 infinite number of estimates of the model which fit the 
data, i.e., which reproduce the data to within the noise 
that is present. Whenever possible, a priori information is 
used to choose which of the estimates to calculate. 
However, a priori information may not be sufficiently 
L5 available or may be suspect. 

In applying the present invention to multiexponential 
decays, an estimate of an unknown relaxation distribution 
(model) is obtained by linearly resolving each point of the 
unknown relaxation distribution as precisely as possible 
2 0 within the limits of the noise. In general, such estimates 
do not reproduce the data. Rather, they are obtained by 
optimizing the linear resolution in a manner that will be 
described later. 

The term "linear" used herein in connection with 
25 "resolution" has the following meaning: 
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Consider a transform operator A which, transforms 
functions (or vectors) xl and x2 to yl and y2 , respectively* 
As equations , this is stated : 

A xl = yl and A x2 = y2 . 
The transform operator A is defined as "linear" if for any 
two real nximbers a and b 

A(a xl + b x2) = a yl + b y2 • 

A transform operator which is a matrix will always have 
linear resolution. It is possible to construct non-matrix 
transform operators which behave similarly to matrices for 
only a restricted set of functions (or vectors) . Moreover, 
a matrix can be derived from a non-matrix transform operator 
that expresses this behavior. 

The present invention involves the following 
relationship between a true model m"^ (y) and data dj^ (e.g., 
multiexponential decays) : 

where I is the appropriate domain of definition and 
gj^(y) represents a data function, more particularly/ one of 
the N data functions which compose the data kernel of a 
linear transform. A data function is a mathematical 
representation of how a model is mapped to a particular data 
point. Equation (1) is referred to as "the forward 
problem. " 



Decaying systems generally have a point, usually 
defined as t = 0 (where t is time) , at which a system begins 
decaying. The system decays to a constant value, often 
zero, as t approaches infinity. In nuclear magnetic 
resonance (NMR) , for example, the t=0 point is when an 
excitation pulse energizes a sample. 

The multiexponential forward problem commonly used in 
data analysis and inverse theory has the form 

In this equation dy, represents the data; m^ represents the 
unknown model, and -e""**/^* represents the data function, 
where r designates the time constant of the exponential 
decay. 

Equation (2) can be expressed in the continuous form 

v^here dQ is the Dirac delta fvinction. 

Incidentally, as used herein, the range of indices such 
as i, j and k are determined by the equation in which they 
are used. For example in the equation 

^iy) = Y:bigi(y) (4) 
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the index i would be assumed to range over all the data 
functions 51(3/) • which are normally numbered with positive 
integers starting at 1. But in the equation 

= (5) 

the index i ranges over the rows of a transform matrix and 
the index j ranges over the columns of the transform matrix, 
which is the same range as the data functions in the 
previous equation. 

It is desirable to plot the relaxation distribution 
(unknown model) versus the logarithm of r, and to determine 
the resolution of a relaxation distribution on a log scale. 
Applying a change of variables y=ln(T) to equation (3), 
without simplification, yields 

/■4-00 
Y, miS(e^ - e»Oe-**'"' {e^dy) ( 6 ) 

Applying the Dirac delta function identity 

to equation (6) along with standard simplification yields 

f^YmAy-yi)e-*'""dy (8) 

J-oo i 

20 Substituting 



L5 
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into equation (8) gives 

dfc= /■^m''(y)e-*'«-|iy- (lo) 

J— OO 

where m^(y) represents the unknovai model, and e""**«~^ 
represents the data function, which may be expressed 

as : 

gk{y)^e'^^^^\ (11) 

Hereinafter, the form of the multiexponential transform in 
equation (10) will be termed the forward problem. 

It should be noted that the data function approaches 1 
as T approaches infinity for all values of tj,. Pursuant to 
the invention, the data function of equation (11) has been 
found to be very useful in multiexponential analysis, but 
other data functions may be appropriate for other 
applications of the invention. 

As stated earlier, a first step in achieving the 
objectives of the invention is to construct a transform 
operator, such as a transform matrix, which maps data (e,g. 
decay signals) to unknown model parameters. Fig. 1 is a 
diagram showing a matrix with coefficients a^^ • y for 
transforming data d^. . . d^ to produce parameters m^ . . . m,, 
of an unknown model. The transform matrix is typically a 
matrix in which each row corresponds to a r of interest. 
Selection of r's of interest and data points for initial 
coefficients will be guided by available information in the 
particular field in which the matrix is to be used. 
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In accordance with the invention, it is desired to 
obtain an estimate of an unknown model with linear 
resolution, and since the data functions used are linear 
functions of the model, the estimate of the unknown model 
can be calculated by multiplying the data by a matrix • The 
matrix is chosen so that each point of the estimate of the 
unknown model linearly resolves the corresponding point of 
the model as well as possible with an acceptable noise gain 
i.e., with optimal linear resolution. Since matrix 
xnultiplication is a linear operation, it yields an estimate 
which does not necessarily reproduce the data, but which 
does have linear resolution. Linear resolution is a 
desirable property of an estimate, because each point of the 
estimate resolves the corresponding point of the model in 
the same way, independent of any particular model. 

A goal in constructing a transform matrix in accordance 
with the invention is to calculate linear combinations of 
the data functions which yield a resolution fimction that 
resolves as small a region of the unknown model as possible. 
Accordingly, it is preferred that each resolution function 
corresponding to a row of the matrix be characterized by 
optimal linear resolution. Preferred criteria for selecting 
coefficients of a transform matrix which yield optimal 
linear resolution are described later. Together, the 
resolution function and its noise gain give a concise 
formulation of the ambiguity of a point in an estimate of an 
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unknov/n model from information given by data and 
corresponding data functions. 

Calculation of a transform matrix proceeds row by row. 
It is convenient to use the same noise gain for all rows of 
a transform matrix, but this is not required. As noted 
above, each row normally corresponds to a r of interest. 
The rows can be ordered by increasing r from top to bottom 
of the matrix. A spacing of 16 r values per decade has been 
found to work well. As an example, r values may be 
calculated between O-lt^^^and iot^^, where t^^^^ is the 
smallest time at which a decay signal is to be sampled and 
^max the largest time. Each row of the matrix corresponds 
to a resolution function that is centered on a t of 
interest . 

It is presently believed that the best way to calculate 
the coefficients of a particular transform matrix (which, 
incidentally, may have only one row) is to solve a 
constrained minimization problem. 

The information needed before the constrained 
minimization can be performed are (1) the data ftinction, 
gi(y)# foT each data point, d^, (2) the standard ^ie via t ion of 
the expected noise, a^, of each data point, (3) the desired 
noise gain, NG, and (4) the r of interest, r^., on which the 
resolution function is to be centered. 

The constrained minimization problem to be solved is to 
minimize I in the equation: 




(12) 
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by varying b^, where 



^(y) = E%.(y) 



(4) 



i 



and where b^ are trial coefficients of a transform matrix 



information on resolution functions can be found in Parker, 
Robert L. 1994; Geophysical Inverse Theory ; Princeton, New 
Jersey; Princeton University Press, incorporated herein by 
reference . 

The trial resolution function preferably complies 
substantially with the following constraints: 

(i) R{y) > 0 [Nonnegative constraint] 

(ii) R(yic) > 1 [Peak constraint] 

(iii) R (y) monotonically decreases from y,^ [Monotonicity 
constraint] 



Constraint (iii) may be satisfied automatically for the 
data function used for the multi exponential problem. Thus, 
while it is stated as a separate constraint, it is to be 
understood that it may, in fact, be satisfied inherently. 



the transform matrix, once the trial coefficients b^ which 
minimize I have been found, it is highly desirable that the 
trial coefficients be scaled by a constant so that the area 



row and R(y) is a trial resolution function. More 




.n constraint] 



To arrive at the final coefficients, a 



for row i of 
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of their associated resolution function is unity on a log 
scale. This is accomplished by 

/+00 ^ 
Y.^k9k{y)dy (13) 

Making the area of the resolution function unity insures 
► that each point in the estimate of the unknovm model is a 
local average. 

Incidentally, for a data point dj, sampled on a decay 
curve at time tj, the data function previously described can 
O be modified by a balancing function B (y) , so that the data 
CO) function becomes 

This modification of the data function has the effect of 
\A making the resulting transform matrix generate an estimate 

■'^5 of the unknown model multiplied by B(y) but with the 

required noise gain. While many balancing functions are 
possible, a useful balancing function if B(y)=^'^ • The 
choice w = 1 may be useful since it increases the amplitudes 
of the relaxation components at late time constants by a 
!0 factor of r . - 
Constructing a transform matrix in which each 
resolution function corresponding to a row of the matrix is 
characterized by optimal linear resolution preferably 
involves a process termed constrained optimization. In the 
2 5 preferred embodiment of the invention, this process includes 
an outer loop, a middle loop and an inner loop. The outer 
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loop converts the constrained optimization problem to a 
series of unconstrained multidimensional minimization 
problems using the simple penalty method described by 
Fletcher (Fletcher, R. 1987 Practical Methods of 
Optimization ; Toronto; John Wiley & Sons) incorporated 
herein by reference. The middle loop converts the 
unconstrained multidimensional minimization problems into a 
series of one dimensional (ID) minimization problems via the 
conjugate gradient method described by Press et al . (Press, 

Teukolsky, S.A., Vetterling, W.T., Flannery, B.P. 
1992; Numerical Recipes in C; The art of scientific 
computing ; 2nd Ed.; Cambridge UK; Cambridge University 
Press) incorporated herein by reference. The 1 dimensional 
minimization problems are solved using the golden section 
search or parabolic interpolation described by Press et al. 

A desirable condition for each resolution function is 
that its area should be s\ibstantially minimized. 
Calculating the area of the trial resolution function I to 
be minimized, given the trial coefficients, bi, can be 
accomplished using the equation: 



where is the area of the data function, gi (y) , given 
by the equation 



(15) 




(16) 
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The integral's upper limit of infinity has been replaced by 
the value Yl, which is a large finite number. Experience 
has shovTn that Yl = In (1000*tj^^) , is a good choice. 

The simple penalty method converts the constrained 
optimization problem to an unconstrained one via residues, 
r^, A residue is a measure of how much a particular 
constraint is violated. If a particular residue is greater 
than zero, then the corresponding constraint is considered 
satisfied. The more negative the residue, the more the 
constraint is violated. The constrained optimization 
problem then becomes the unconstrained one of minimizing 
in the equation : 

■^'' = i:*i-^ + ^i:"»in(rm,0)* (17) 

ft m 

where min(x,y) returns the minimxim and x and y. The 
value of P starts small, and is then increased by a small 
factor, for example, 0.1, after each unconstrained 
minimization. The solution to one constrained minimization 
is used as a starting point for the next constrained 
minimization. The value of P is increased until has 
stabilized. The test for stability is whether the value of 

changes by an accumulative factor of 10"^ over four 
consecutive increases in P. 

Each of the non-negative, peak, monotonicity and noise 
gain constraints will usually have associated residues. 
Before the residues can be calculated, the data functions 
and the resolution functions are discretized. 
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The data functions and the resolution functions are 
discretized by sampling them at e.g., 16 points per decade 
of r starting at t= 0.01 time units and continuing to t= 
1000*tjj^3, time units. The discretized data functions, gi (y) , 
are represented as g^^ . The total niomber of points at which 
the functions are sampled is denoted L. The discretized 
trial resolution function is represented by R^. 

Thus 

The residues 
defined to be 

giving L residues. 

The residues of the peak constraint, r^, are defined to 

be 

where Ip is the index of the discretized trial 
resolution function at the r of interest. This yields one 
residue . 

The residues of the monotonicity constraint are defined 

to be 

^MJ^^^„Ki-i 2>l>lp (21) 

and 



Rt^'£^bi9u (18) 



of the nonnegative constraint, t-^^ , are 



rr = JRt (19) 



(22) 



These two equations together yield another L - 1 residues. 

The residue of the noise gain constraint is defined to 

be 



This equation yields one additional residue. 

In total there are 2*L+1 residues to express the four 
constraints . 

The conjugate gradient technique is explained by Press 
et al . The termination conditions used may be an 
accumulative factor change in of no more than 10"^ over 25 
consecutive iterations. The total number of iterations can 
be limited to 3000*N, where N is the niimber of data points. 

The one dimensional minimization problems are solved 
using the golden section search or parabolic interpolation 
described by Press et al • The ID optimization is terminated 
if the functional value does not decrease by a certain 
fraction each step. The fraction can be set to 10"®, for 
example. A step limit of 1000 on each ID optimisation can 
be set to prevent infinite or unproductive near infinite 
loops . 

The linear transform for the trial row coefficients to 
the trial resolution function would usually be computed a 
large number of times during each one dimension optimization 
and would be computationally demanding. However, since the 
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relationship between the trial row coefficients and the 
trial resolution function is always linear, this linearity- 
can be used to greatly reduce the ntimber of times the trial 
resolution functions have to be calculated directly from the 
5 trial coefficients. 

In one dimensional optimization, a scalar parameter c 
is varied along a line with direction d^; the direction 
being provided by the conjugate gradient method. Therefore, 
for the trial coefficients b^, optimization occurs along a 
^^0 line defined by 

bn^b^+cdf. (24) 

where b° is also supplied by the conjugate gradient method. 

Because of the linearity relationship between b^ and R^, a 
d similar equation can be written for the corresponding trial 

^riS resolution function R^^, 

[i Bi^mi + cdf (25) 

where 

^^^a9gu (26) 
2 0 * 

and 

'^f^JZdtga. (27) 

i 

Therefore, along each line of optimization, the linear 
transform from the trial coefficients to the trial 
25 resolution function only needs to be calculated twice, to 
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calculate R-^ and d^, instead of once for every value of c 
considered. This results in a many fold reduction in the 
one dimensional optimization time . 

The effectiveness of the minimization can be checked by 
how close the value of the peak of the resolution function 
is to unity, since the minimvim area should yield a peak 
value of exactly unity. 

An upper limit on the noise gain, NG, of a particular 
point in an estimate can be applied by the noise gain 
constraint stated earlier, where the area of the resolution 
function is included because the area of the resolution 
function must be normalized (set equal to 1) before the 
noise gain is calculated. 

Figure 2 is a flow chart for calculating the 
coefficients of a row of a transform matrix pursuant to the 
foregoing description, using a digital computer 
conventionally. Coefficients for successive rows of the 
matrix can be calculated in the same manner for each t of 
interest and for each noise gain deemed appropriate. For 
calculating the coefficients of any row, the coefficients of 
the preceding row can be used as starting points... 

Constrained optimization methods may perform better if 
all the data have noise with a standard deviation of 1. 
Fortunately, it is quite simple to modify the data functions 
so this is the case. Given the real data functions, gi (y ) # 
the adjusted data function, gi*(y) is 

i^(y) = ffi(i/)M (28) 



The adjusted data functions are then supplied to the 
constrained optimization method along with the assumption 
that the standard deviations of noise at all the data points 
are 1. The final coefficients, a^^, produced by the 
constrained optimization based on the adjusted data 
functions, need to be corrected for the adjustment. The 
correction to the final coefficients is 

<Hj==<^ijM (29) 

If a data set has 100' s, 1000' s, or even more evenly 
spaced points on a decay curve, it is much more efficient 
computationally for calculating the coefficients of the 
transform matrix (and also for applying it to data) to 
average adjacent data points together to create a new group 
data point. The corresponding data functions must also be 
averaged together to get the grouped data function 
corresponding to the new grouped data points . The size of 
the groups is important* While, in general, it is best to 
have larger groups at later sample times, groups which are 
too large will reduce the resolution of the resolution 
functions. Groups which are too small are inefficient. A 
logarithmic group size appears to be a good choice. 
The equation 

Gn = maxihint{C 10^'^). (30) 

appears to produce good group sizes for C = 0.1, D = 10, 
where 1 is a positive integer starting at 1 and increasing 



to whatever value is appropriate. The function max(x,y) 
returns the maximuin of x and y, and the function int(2) 
rounds z down to the nearest integer* For 512 data points 
the above equation gives group sizes of 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 3 3 5 6 7 10 12 15 19 25 31 39 50 63 79 100 30 

Occasionally, instead of measuring data at a point in 
time on a decay curve, data will be measured by averaging 
over a window between two points in time. In this case, the 
data function can be approximated by averaging together a 
larger nximber of sample points over the window. 

Modification of the data functions will also be 
required if a trigger that initiates a decay curve is not a 
single point in time. For example, in time resolved 
fluorescence spectroscopy, the flash of light triggering the 
decay will last a finite length of time. If the intensity 
versus time for the flash is L{t) then each data point will 
be convolved with this function. The corresponding data 
functions will also have to be convolved with the same 
function. 

Fig. 3 illustrates a transform matrix constructed in 
accordance with the invention. Three sets of coefficients 
for T=l, 10 and 100 are shown as columns so that they fit on 
the page. Each set has 48 coefficients calculated from 48 
data functions. The first 32 data functions correspond to 
32 sampling points that are evenly spaced by one time unit 
at times 1 to 32. The next 16 data functions are evenly 
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spaced by 30 time units between time 62 and time 512. Tlie 
standard deviation of the noise for the first 32 points is 
assumed to be 1, and the standard deviation of the noise of 
the next 16 is assumed to be 0.182 57. The standard 
> deviation of the noise of the last 16 is a factor of 
l/sqrt(3 0) less than the first 32. This drop in the 
standard deviation of the noise would result if the cut-off 
frequency of a low pass filter before the analog- to-digital 
converter were dropped by a factor of 3 0 before the point 
'^:b was acquired- For the sake of simplicity, the matrix was 
S constructed without the use of balancing functions or data 

=□ function groups. 

y Fig. 4 illustrates data functions, resolution 

functions, noise response, and point spread functions for a 
IjiS matrix constructed in accordance with the invention where 
^ the noise gain is 1.000. The data function curves are for 

time samples at t=l, t=2 . . . t=N from left to right. Each 
resolution function corresponds to a set of data functions. 
The abscissa in each diagram is in xinits of t on a log scale 
20 [ln(T)], and each resolution function is localized about a 

particular r value. As a result, each point spread, function 
tends to be localized about a particular t value. 

Fig. 5 is a diagram similar to Fig. 4 but for a matrix 
constructed with a noise gain of 3.1623. It should be noted 
2 5 that in each of Figs. 3, 4, and 5 there are 4 8 data points 
(time samples) . 
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A comparison of resolution functions produced by 
matrices with different noise gains is shown in Fig. 6. 
Higher resolution is achieved with higher noise gains, but 
there is a trade-off between resolution and noise gain. 
Greater noise tends to make the results achieved less 
reliable . 

In addition to resolution functions, performance of a 
transform matrix can be judged using the PSF' s and noise 
response. To judge the performance of a transform matrix, 
information is required as to corresponding time points on 
the decay curve at which data should be acquired as well as 
assumed noise at each data point. 

Fig. 7 shows point spread functions associated with 
four matrices constructed in accordance with the invention 
for different noise gains. PSF's can be calculated using 
simulated data at required time points for monoexponential 
decays at r = 1,2/5,10,20,50,100, ... time units, assuming 
the initial amplitude of the decay curve is 1. No noise is 
added to the simulated data. It should be noted that in 
each of Figs. 6 and 7, and also on Figs. 9 and 10 to be 
described, there are only 32 equally spaced data-points 
(with the same noise) . 

To calculate the noise response, several realizations 
of a decay curve which are purely noise are generated. The 
noise may be assumed to be Gaussian. The standard deviation 
of the generated noise is scaled so that the standard 
deviation of the noise of the first data point is 1. Then 
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the sampled noise decay curves are multiplied by th.e 
transform matrix to obtain the relaxation distribution. 
Several realizations of the relaxation distribution of the 
noise can be plotted to obtain a "feel" for the 
distribution. 

An important characteristic at each point in the 
estimate of the unknovzn model is the standard deviation due 
to the noise in the data. If the noise in the data is 
uncorrelated, has a mean of zero, and all points have the 
same finite standard deviation, it can then be characterized 
by a single standard deviation. The noise gain for each 
point can be defined to be the standard deviation of the 
point divided by the standard deviation of the noise in the 
data and can be calculated directly from the coefficients 
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Once a transforro matrix has been constructed, 
experiments can be performed in which the spacing of the 
data acquisition points, the number of data acquisition 
points, and the grouping of data acquisition points are 
changed* The results of these experiments can be cpmpared 
by observing the resolution functions and/or the PSFs that 
are produced. As a rule of thumb, the linear resolution is 
proportional to the height of the resolution function 
provided that the resolution function has \inity area on a 
log scale. 



Each trial resolution function is localized about a 
point of interest in the unknown model by requiring the peak 
of the resolution function to have a value greater than 
unity at a r of interest and then by minimizing the area of 
the resolution function. The effectiveness of the 
minimization can be checked by how close the value of the 
peak of the resolution function is to unity, since the 
minimxim area should yield a peak value of unity. 

From the foregoing description, it is evident that each 
transfoorm matrix is constructed so that each point of the 
estimate of the unknown model linearly resolves the 
corresponding point of the unknown model as well as possible 
within an acceptable noise gain. Since matrix 
multiplication is a linear operation, it yields an estimate 
which does not necessarily reproduce the data but does have 
linear resolution. With linear resolution, each point of 
the estimate resolves the corresponding point of the unknown 
model in the same way independent of any particular unknown 
model, A highly desirable property of linear resolution in 
providing an estimate of an unknown model is that it enables 
a human interpreter to obtain an intuitive "feel/' of what 
the data reveals and does not reveal about the unknown 
model . 

Each resolution function gives a concise mathematical 
description of the linear resolution at each point of the 
unknown model and is independent of the unknown model and 
the data. Viewing the transform matrix as a digital lens, 
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linear combinations of the data functions are calculated 
which yield a resolution function that resolves as small a 
region of the unknown model as possible. 

After a transform matrix is constructed and selected, 
it is incorporated in a signal processor of a computer, as 
software or hardware, for example, as indicated in Fig, 8* 
In this figure the INPUT represents a source of sampled 
digital signals such as multiexponential decays, e.g., NMR 
decays obtained from well - logging , The DATA PROCESSOR and 
STORAGE are components of a conventional digital computer. 
The OUTPUT MODULE may have conventional DISPLAY, NETWORK 
INTERFACE, and PRINTER COMPONENTS, for example. 

An estimate of an unknown model is calculated by 
multiplying data, e.g., multiexponential decays, by the 
transform matrix. Several transform matrices for different 
noise gains, for example, can be incorporated in the signal 
processor and accessed selectively to provide a user with 
greater flexibility. 

As stated earlier, an object of the present invention 
is to provided a better estimate of an unknown model, from 
which useful information can be obtained. In the 
aforementioned Prammer patent. Figure 8 of the patent is a 
mapping of estimated NMR signal decay times into pore sizes 
of an investigated earth formation. The curve of Fig. 8 is 
an estimated relaxation distribution. The present invention 
provides a better estimate of the relaxation distribution, 
and also a better signal-to-noise ratio (SNR) . Resolution 
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functions produced by matrices constructed in accordance 
with the invention are superior, in terms of peak value and 
localization to resolution functions produced by matrices of 
the Prammer patent for the same data and noise. 

In general, the present invention can be used to 
provide better estimates of an unknovm model than the prior 
art, estimates from which more reliable and useful 
information can be obtained. These improved results are 
achieved by emphasizing linear resolution irrespective of 
whether data fits a model, and, in fact, without any attempt 
to fit data to a model. The invention is particularly 
useful in multiexponential signal processing, such as in the 
analysis of multiexponential decays. 

As mentioned earlier, with particular reference to the 
REVIEW article, multiexponential analysis is useful in a 
host of scientific and technological applications. One of 
those applications, NMR, such as NMR in well -logging, has 
already been discussed. Another application is MRI 
(magnetic resonance imaging) . Fig. 9 shows decay curves in 
an MRI application. The decay curves represent pixels taken 
from a series of 32 magnetic resonance images of., the brain 
of a multiple sclerosis patient. The T2 relaxation data 
were acquired with a 10ms sample time out to 320ms. The 
strength and decay of the signals contain valuable 
information about the tissues. Some of the major tissue 
types of interest in multiple sclerosis are normal appearing 
white matter (NAWM) , cerebral spinal fluid (CSF) and 



lesions, which can be classified as chronic or acute. Fig, 
10 shows relaxation distributions yielded by applying to the 
decay curves of Fig. 9 transform matrices of the invention 
with various noise gains. The larger the noise gain of the 
estimate, the wider the range of relaxation estimates 
available. This is a characteristic which shows up in the 
resolution functions and is due to the fact that the larger 
the noise gain the more likely there is a feasible 
resolution function available. 

Estimating the noise in the estimates shown in Fig, 10 
can be accomplished in several ways. The first is to 
estimate the noise in the data and multiply by the noise 
gain for each transform matrix. The ideal way to measure 
the noise in the data is to repeat the measurements a large 
number of times and calculate mean, standard deviation and 
covariance. These statistics can then be propagated through 
transform matrices using standard statistical procedures. 
Unfortunately, the measurement of the decay curves takes 
about 2 0 minutes to complete on a patient, so large numbers 
of repetitions are impractical. 

Another way to measure the noise in the data is to 
estimate the standard deviation from previous measurements 
using the same instrxoment. This is not always reliable 
because the noise can vary from sample to sample. For the 
data in Fig. 9, previous measurements of noise gives the 
standard deviation of about 10 scanner units. After 
multiplying the standard deviation by the noise gain, the 



predicted standard deviation for the noise in Fig. 10 is 
100. 

A third way to estimate the noise is to consider that 
while the noise gain increases by a factor of 10 in Fig, 10, 
5 the resolution gain increases by only 1.4 for relaxation 
rates around five sample times. In Fig. 10, a portion of 
the signal between 3 0 and 2 0 0ms increases proportionately to 
the noise gain. This strongly suggests that the portion is 
due to random uncorrelated additive noises in the data. It 
izD is possible then to measure the standard deviation of the 
m noise between 3 0 and 2 0 0ms and work back to the noise in the 
:.g data. 

Further applications of the invention will be apparent 
7' to persons of ordinary skill in the art. For example, 

:';5 decaying sinusoids are a common problem in inverse theory. 

A decaying sinusoid can be handled if it is band pass 
3 filtered at the particular bandwidth of interest and then 

the magnitude of the decay curve taken. The relaxation 
distribution of the magnitude decay curve can then be 
20 calculated. 

Any row of a transform matrix, since it corresponds to 
a resolution function, can be applied to a time series in 
the same way as digital filters. The "relaxation" digital 
filter will be useful for applications such as ultrasound 
2 5 and radar. Using quadrature detection, it would be possible 
to measure the magnitude of a reflection off of an 
interface. If the signal from an interface of interest 
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oscillated for a while after the initial sound wave had 
passed, the decay time of the oscillation would reveal 
information about the interface. Applying various 
relaxation digital filters to the ultrasound time series 
) would allow characterization of the interfaces. If it were 
desirable to detect a particular range of relaxation time, 
selection of the coefficients of a relaxation digital filter 
would give a resolution function with a more boxcar- like 
shape . 

4p Other applications of the invention include 

photoliaminescence and time -resolved fluorescence 
=D spectroscopy of biological and other types of samples, 
H electrical signals radiating from ore bodies in geophysical 

exploration and acoustic, electrical and electromagnetic 
'Zj5 decay processes. It is possible to design a transform 
^ni matrix which handles data integrated over a window or 

unevenly spaced in time by modifying the data functions. 

Principles of the invention can be applied to problems 
which data functions other than decay curves. If the data 
20 functions are cosine functions, resolution functions can be 
generated for low pass, band pass and high pass filters as 
well as windows for discrete Fourier transforms. A low pass 
filter can be designed by requiring maximxim area below the 
cutoff frequency, minimum area above the cutoff and an 
25 optional requirement of monotonicity to eliminate wiggles. 

The bounds on all parts of the filter would be 0 and 1. A 
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limit on broadband noise gain could also be imposed to 
improve the robustness of the filter. 

While preferred forms of the invention have been shown 
and described, these forms are intended to be exemplary, not 
5 restrictive, and it will be apparent to those skilled in the 
art that modifications can be made without departing from 
the principles and spirit of the invention, the scope of 
which is set forth in the appended claims. 
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WHAT IS CLAIMED IS: 

L 1. A computer- readable mediuia containing a transform 

I operator constructed to provide an estimate of an unknovm 

3 model with substantially optimal linear resolution. 

1 2- A computer-readable medium according to Claim 1, 

2 wherein the transform operator comprises a matrix having at 

3 least one row of coefficients corresponding to a resolution 

4 function. 

3, A computer- readable medixim according to Claim 2, 

S wherein the resolution function has sxibstantially no 

3 negative values . 

%ll 4. A computer- readable medium according to Claim 2, 

iJ2 wherein the resolution function has a peak value of at least 

U4:3 substantially 1 substantially at the point of interest. 

1 5. A computer- readable medium according to Claim 2, 

2 wherein the resolution function monotonically decreases from 

3 its peak value. 

1 6. A computer- readable medixim according to Claim 2, 

2 wherein the resolution function has substantially no 

3 negative values, has a peak value of at least substantially 

4 1 substantially at the point of interest, and decreases 

5 substantially monotonically from its peak value. 
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7. A computer-readable medium according to Claim 6, 
wherein the resolution function complies substantially with 
the following noise gain constraint: 



8. A computer -readable medium according to Claim 2, 
wherein the matrix has a plurality of rows of coefficients, 
each row corresponding to a different resolution function. 

9. A computer -readable medium according to Claim 8, 
wherein the matrix is constructed for use in 
multiexponential decay signal processing and wherein each 
resolution function is substantially centered on a different 
time constant. 

10. A computer -readable medium according to Claim 9, 
wherein each coefficient has a corresponding data function 
in the form c"****"*' or derivations of this form, where y is 
the natural logarithm of a time constant and t^. is a 
sampling point on the decay signal, 

11. A computer- readable medi\im according to Claim 2, 
wherein the coefficients are selected so that the 
corresponding resolution function has a peak value of at 
least substantially 1 substantially at the point of interest 
and wherein the area of the resolution function is 
substantially minimized. 




[IQolse gain constraint] 



12, A computer- readable mediiam according to Claim 11, 

I wherein the area of the resolution function is set to be 

I substantially 1 on a logarithmic scale. 

1 13 . A method of constructing a transform operator in 

2 which a plurality of coefficients are selected to produce a 

3 corresponding resolution function that provides 

4 substantially optimal linear resolution. 

■^1 14. A method according to Claim 13, wherein the 

-""2 resolution function has sxibstantially no negative values, 

.113 has a peak value of at least substantially 1 substantially 

Q 4 at the point of interest, decreases substantially 

; - 5 monotonically from the peak value, and has a gain within a 

: J 6 predetermined range . 

1 15. A method according to Claim 13, wherein the area 

2 of the resolution function is set to be siibstantially 1 on a 

3 logarithmic scale. 

1 16 . A method according to Claim 13, wherein the 

2 transform operator is constructed to provide substantially 

3 optimal linear resolution between data and an unknown model 

4 without any consideration of whether the data fits the 

5 unknown model • 
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17. A method of multiexponential signal processing, 
which comprises: 



1 sampling a multiexponential signal and applying the 

t transform operator of Claim 1 to the sampled signal. 

1 18. Apparatus for multiexponential signal processing, 

2 which comprises a signal processor that has the transform 

3 operator of Claim 1. 

1 19. A method of constructing a transform operator, in 

2 which a plurality of coefficients are calculated using data 

3 functions that emphasize linear resolution irrespective of 

4 data. 

1 2 0. A method according to Claim 19, wherein each data 

2 function is in the form e"*"**"* or derivations of this form, 

3 where y is the natural logarithm of a time constant and tj, 

4 is a sampling point on a decay signal, 

1 21. A method of exponential signal processing, which 

2 comprises: 

3 providing a sampled multiexponential signal; and 

4 applying the sampled multiexponential signal to a 

5 transform operator constructed to provide substantially 

6 optimal linear resolution between outputs of the transform 

7 operator and an unlmown model . 
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ABSTRACT OF THE DISCLOSURE 

For signal processing of transients such as 
mul tiexponential decays, a transform operator (e.g., matrix) 
is constructed by emphasizing linear resolution, rather than 
using fitting routines that attempt to find an estimate of 
an unknown model that fits data* Use of the transform 
operator to process multiexponential signals produces 
outputs that are a better estimate of the unknovTn model or 
of some segment of the unknown model* 
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FIG. 7 
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FIG. 8 
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FIG- 10 
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